##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(ggsci)
library(RColorBrewer)
library(ggpubr)
library("scales")

##########################################################################################

option_list <- list(
    make_option(c("--info_file"), type = "character") ,
    make_option(c("--public_file"), type = "character") ,
    make_option(c("--images_path"), type = "character")
)

if(1!=1){

    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    info_file <- paste(work_dir,"/baseTable/STAD_Info.addBurden.MSI_MSS.addCNVType.tsv",sep="")
    public_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/combine/MutationInfo.combine.addMolecularSubType.rmMIX.tsv"
    images_path <- paste(work_dir,"/images/mutBurden",sep="")
}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

public_file <- opt$public_file
info_file <- opt$info_file
images_path <- opt$images_path

dir.create( images_path , recursive = T )

###########################################################################################

dat_sample <- data.frame(fread( info_file ))
dat_public <- data.frame(fread( public_file ))

###########################################################################################
## NJMU的多个样本负荷取中位数
dat_tmp_nmu <- dat_sample %>%
group_by( Patient , Age , Class , Tobacco , Alcohol , HP ) %>%
summarize( BurdenExon = median(BurdenExon) )
dat_tmp_nmu <- data.frame(subset(dat_tmp_nmu , Class == "IM"))

## 公共数据的样本
dat_tmp_public <- dat_public[,c("Tumor" , "Age" , "Class" , "Tobacco" , "Alcohol" , "HP" , "BurdenExon" , "Molecular.subtype")]
colnames(dat_tmp_public)[1] <- "Patient"

## 注释分子亚型
dat_tmp_nmu <- merge( dat_tmp_nmu , dat_tmp_public[,c("Patient" , "Molecular.subtype")] , by = "Patient" )

## 使用样本
dat_plot <- rbind( dat_tmp_nmu , dat_tmp_public )
dat_plot$Tobacco <- factor( dat_plot$Tobacco , levels = c("No" , "Smoke") )
dat_plot$Alcohol <- factor( dat_plot$Alcohol , levels = c( "No" , "Drink") )
dat_plot$HP <- factor( dat_plot$HP , levels = c("Negative" , "Positive") )

## 提取年龄信息全的样本
dat_plot <- subset( dat_plot , !Age %in% c("[Not Available]" , "unknown") )
dat_plot$Age <- as.numeric(dat_plot$Age)
dat_plot$Molecular.subtype <- ifelse( dat_plot$Molecular.subtype %in% c("MSI" , "POLE") , "MSI/POLE" , dat_plot$Molecular.subtype )
dat_plot$Molecular.subtype <- ifelse( dat_plot$Molecular.subtype %in% c("EBV" , "unknown") , "unknown" , dat_plot$Molecular.subtype )

if(1!=1){
    ## 制造MSS
    dat_plot2 <- subset(dat_plot , Molecular.subtype != "MSI/POLE")
    dat_plot2$Molecular.subtype <- "MSS"
    dat_plot <- rbind( dat_plot , dat_plot2 )
}

###########################################################################################
## 按照分子亚型对IGC和DGC进行比较，矫正年龄
#mol_list <- c( "MSS" , "GS" , "CIN" , "MSI/POLE")
mol_list <- c( "GS" , "CIN" , "MSI/POLE")

result <- c()
for( molN in mol_list ){
    # logistic回归模型拟合,矫正年龄
    print(molN)
    dat_plot$use_col <- dat_plot[["Molecular.subtype"]]
    dat_plot$useCol <- dat_plot[["BurdenExon"]]
    
    tmp_dat <- subset( dat_plot , Class %in% c("IGC" , "DGC") & use_col == molN )
    tmp_dat <- tmp_dat[!is.na(tmp_dat$use_col),]
    model <- glm( useCol ~ Class + Age , data = tmp_dat)

    ## p值
    P <- summary(model)$coef[2,4]
    ## 95%可惜区间
    OR <- cbind(OR = coef(model), confint(model))[2,]
    tmp <- data.frame( baseType = molN , 
        model = "adjust_age" , 
        var = "GC" , 
        time = "IGC vs DGC" , p = P , coef = OR[1] , coef_025 = OR[2] , coef_975 = OR[3] )

    result <- rbind(result , tmp)
}

result_gc <- result

###########################################################################################
## IM中比较饮酒，分不同分子亚型
base_use <- "Alcohol"
result <- c()
## 单因素
# logistic回归模型拟合
print(base_use)
dat_plot$use_col <- dat_plot[[base_use]]
dat_plot$useCol <- dat_plot[["BurdenExon"]]

## 分不同分子计算
for( molN in mol_list ){
    tmp_dat <- subset( dat_plot , Class %in% c("IM") & Molecular.subtype == molN )
    tmp_dat <- tmp_dat[!is.na(tmp_dat$use_col),]
    model <- glm( useCol ~ use_col + Age , data = tmp_dat)

    ## p值
    P <- summary(model)$coef[2,4]
    ## 95%可惜区间
    OR <- cbind(OR = coef(model), confint(model))[2,]
    tmp <- data.frame( baseType = molN , 
        model = "adjust_age" , 
        var = "Alcohol drinking" , 
        time = "Yes vs No" , p = P , coef = OR[1] , coef_025 = OR[2] , coef_975 = OR[3] )
    result <- rbind(result , tmp)
}

result_im <- result
result <- rbind(result_gc , result_im)

out_name <- paste0( images_path , "/mutBurden_logistic_adjustage.tsv" )
write.table( result , out_name , row.names = F , quote = F , sep = "\t" )
